using System;

using Atomic.Libraries;

namespace Atomic.Thermodynamics.Phonons
{
	[Serializable]
	public class PhononPotential : IPotential
	{
		[NonSerialized] private RollingCache<double, Tuple<double, double, double>> values;

		public PhononPotential(IPhononStateDensity stateDensity, int atoms)
		{
			if (stateDensity == null)
			{
				throw new ArgumentNullException();
			}

			if (atoms < 1)
			{
				throw new ArgumentOutOfRangeException();
			}

			StateDensity = stateDensity;
			Atoms = atoms;
		}

		public double FreeEnergy(double temperature)
		{
			// Currently per primitive cell!!!
			return Compute(temperature).Item1;
		}

		public double Entropy(double temperature)
		{
			return Compute(temperature).Item2;
		}

		public double HeatCapacity(double temperature)
		{
			return Compute(temperature).Item3;
		}

		private Tuple<double, double, double> Compute(double t)
		{
			if (values == null)
			{
				values = new RollingCache<double, Tuple<double, double, double>>(10000);
			}

			Tuple<double, double, double> value;
			if (!values.TryGetValue(t, out value))
			{
				double f, s, cv;
				Integrate(t, out f, out cv, out s);

				// Cache computed free energy, entropy, and heat capacity as a tuple.
				value = Tuple.Create<double, double, double>(f, s, cv);
				values.Add(t, value);
			}

			return value;
		}

		private void Integrate(double t, out double f, out double cv, out double s)
		{
			double hbar = 6.5821192815e-16; // eV*s
			double kb = 8.617332478e-5; // eV/K

			double d = 0.0;
			f = 0.0;
			cv = 0.0;
			s = 0.0;

			double omega0 = 0.0;
			double di0 = 0.0;
			double fi0 = 0.0;
			double cvi0 = 0.0;
			double si0 = 0.0;
			foreach (IPhononStateDensityPoint point in StateDensity.Points)
			{
				// Change to SI units; convert frequency to angular frequency.
				double omega = 2.0 * Math.PI * point.Frequency;
				double g = point.Density / (2.0 * Math.PI);

				double a = hbar * omega / (2.0 * kb * t);
				double a1 = Math.Sinh(a);
				double a2 = a / a1;
				double a3 = Math.Log(2.0 * a1);

				double di = g;
				double fi = g * a3;
				double cvi = g * a2 * a2;
				double si = g * (a / Math.Tanh(a) - a3);

				if (omega0 > 0.0)
				{
					// Ignore non-positive (negative, imaginary, ...) frequencies. Need to check left-most end-point only.
					// More importantly, this also implicitly skips the first point for trapezoid evaluation.
					// Approximate integrals by trapezoids.
					double b = 0.5 * (omega - omega0);
					d += b * (di + di0);
					f += b * (fi + fi0);
					cv += b * (cvi + cvi0);
					s += b * (si + si0);
				}

				// Use end-points for next trapezoids.
				omega0 = omega;
				di0 = di;
				fi0 = fi;
				cvi0 = cvi;
				si0 = si;
			}

			// The density of state is per primitive cell. The Brillouin zone is also with respect to the primitive cell.
			// Thus these are also per primitive cell.
			f *= kb * t; // eV
			cv *= kb; // eV/K
			s *= kb; // eV/K

			// Scale to the specified structure size. Three phonon states per atom.
			double c = 3.0 * Atoms / StateDensity.PhononStates;

			//if (normalize)
			//{
			//	// Correction factor.
			//	c *= StateDensity.PhononStates / d;
			//}

			f *= c;
			cv *= c;
			s *= c;
		}

		public IPhononStateDensity StateDensity
		{
			get;
			private set;
		}

		/// <summary>
		/// The number of atoms in the structure of interest. This is used for scaling of energies if the size of
		/// the structure is different from the size of the primitive structure.
		/// </summary>
		public int Atoms
		{
			get;
			private set;
		}
	}
}
